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Abstract 

We study the zero-temperature equation of state (EOS) of solid 4 He in the hexagonal closed 
packet (hep) phase over the — 57 GPa pressure range by means of the Diffusion Monte Carlo 
(DMC) method and the semi-empirical Aziz pair potential HFD-B(HE). In the low pressure regime 
(P ~ 0—1 GPa) we assess excellent agreement with experiments and we give an accurate description 
of the atomic kinetic energy, Lindemann ratio and Debye temperature over a wide range of molar 
volumes (22 — 6 cm 3 /mol). However, on moving to higher pressures our calculated P — V curve 
presents an increasingly steeper slope which ultimately provides differences within ~ 40% with 
respect to measurements. In order to account for many-body interactions arising in the crystal 
with compression which are not reproduced by our model, we perform additional electronic density- 
functional theory (DFT) calculations for correcting the computed DMC energies in a perturbative 
way. We explore both generalized gradient and local density approximations (GGA and LDA, 
respectively) for the electronic exchange-correlation potential. By proceeding in this manner, we 
show that discrepancies with respect to high pressure data are reduced to 5 — 10% with few 
computational extra cost. Further comparison between our calculated EOSs and ab initio curves 
deduced for the perfect crystal and corrected for the zero-point motion of the atoms enforces the 
reliability of our approach. 



PACS numbers: 67.80.-s,61.50.Ah,61.50.Ks 



I. INTRODUCTION 



The physics of helium at low temperatures is among the most intriguing and intensively 
studied topics in condensed matter science. Despite of being a rare gas element with one of 
the simplest possible electronic structures, helium constitutes a fundamental system which 
is challenging for the test and development of methods based on quantum theory. Because 
of its light atomic mass and weak interatomic interaction, helium is the only system that 
remains liquid under its own vapor pressure (P = 0) at zero temperature. Below 2.17 K, 
liquid 4 He features superfluidity and Bose-Einstein condensation, two striking and inherent 
quantum effects. With an external pressure of ~ 25 bar, the fluid at T = crystallizes into 
the hexagonal closed packet structure (hep), which remains the stable phase of solid helium 
at T ^ and high pressures, made the exception of an fee loop along melting in between 
15 - 285 KMS 

Solid helium is manifestly a quantum crystal. In the regime of ultralow temperatures 
(few mK) this system possesses extraordinarily large atomic kinetic energy (E^ ~ 24 K) 
and Lindemann ratio (7 ~ 0.26), and likewise anharmonic effects on it are of relevance for 
predicting and understanding its thermodynamic properties. 4 Further testimony about the 
uniqueness of this solid is posed by the long-standing controversy sparked by recent experi- 
mental findings about whether perfect crystalline 4 He may exhibit superfluid-like behavior 
and Bose-Einstein condensation (supersolid).-^^ From a technological side, solid helium 
also has some relevance since it is considered as the best quasi-hydrostatic medium hence 
modern technologies based on it have emerged and induced considerable progress in the field 
of high-pressure experiments.— iiLi2ii2 

In this paper, we study the zero-temperature equation of state of bulk solid 4 He in 
the hep phase over a wide pressure range (0 — 57 GPa) with the Diffusion Monte Carlo 
(DMC) method and the HFD-B(HE) Aziz pair potential (hereafter referred to AzizII),— 
and additionally with electronic density functional theory (DFT) to account for many-body 
interactions arising in the system with increasing pressure. In all the work we differentiate 
between two pressure regimes, namely low pressure (0 < P < 1 GPa) and high pressure (1 < 
P < 57 GPa). Quantum Monte Carlo (QMC) methods have proved among the most accurate 
and reliable tools for solving quantum many-body problems associated to condensed matter 
systems.— ii&iLi& In particular, the DMC method is a zero-temperature approach which 
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yields exact estimation (only subject to statistical uncertainty) of the ground-state energy 
and related properties of many-boson interacting systems.— i22i2I During the last few decades 
this and other Monte Carlo techniques (mainly, the variational Monte Carlo -VMC- and Path 
Integral Monte Carlo -PIMC- methods) have been fruitfully applied to the study of noble 
gases and light elements and compounds like He, Ne, H, D, LiH and LiD in homogeneous 
and inhomogeneous phases and both in bulk and in reduced dimensionalities.— i22iMi2£i2£i2L2& 
The great capability of the DMC method is related to the existence of accurate interatomic 
potentials, which are expressed in the form of many-body expansions, and are tuned to 
reproduce empirical and/or theoretical data. Interatomic potentials are of precious value 
because provide computational affordability by allowing one to model atoms as interacting 
points (thus avoiding direct treatment of the electronic degrees of freedom of the system), 
and also simplified understanding of the system under study. In the case of helium, the 
semi empirical pair-potential HFD-HE2 proposed by Aziz et al— more than twenty years 
ago has allowed for quite precise reproduction of the energetic and structural properties of 
the liquid and solid phases near equilibrium.— 1 ^ In this work, we use a newer version of this 
potential, namely the HFD-B(HE) one, 14 which has demonstrated excellent performance in 
the description of the liquid 3 ^ but heretofore has not been tested in the crystal upon high 
pressure. 

Anticipating some of the outcomes we are to present, excellent agreement between our 
results for EOS and experiments is assessed in the low pressure regime for volumes ranging 
from V = 21.30 cm 3 /mol to V = 8.50 cm 3 /mol; however, discrepancies start to develop 
at smaller volumes (P > 0.65 GPa). Within the low pressure regime, we provide exact 
estimation within some statistical error of the kinetic energy per atom, Lindemann ratio and 
Debye temperature of the system by means of the pure estimator technique.— In the high 
pressure regime, however, our equation of state systematic and increasingly overestimates 
the pressure. Discrepancies with respect to measurements amount to ~ 10% at P = 1 GPa 
and to ~ 40% at P = 57 GPa. Previous PIMC work on the EOS of solid 4 He at ambient 
temperature (T = 300 K) and performed with akin model pair potentials arrived to similar 
disagreements.— 1 ^ With the aim of analyzing the possible causes of this large disagreement 
we first examine the influence of finite size effects in our results. Indeed, finite size effects 
become larger by increasing pressure because cut-off distances involved in the calculation 
of the atomic interactions within the system are continously reduced (generally these are 
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chosen as half the length of the simulation box). Accordingly, the radial pair distribution 
function for crystals, g(r), emerges progressively less smooth with compression. Therefore, 
customary corrections devised for dealing with finite size effects which are based on simple 
approximations for g(r), might introduce appreciable deviations in the results (see Fig. [[]). 
Because of these effects, we have regarded as essential to quote the energy tails accounting 
for the finite size corrections by means of two different approaches: (i) considering g(r) ~ 1 
beyond a certain cut-off distance and then integrating the simplified analytic expressions for 
the tails, and (ii), computing the variational Monte Carlo energy of progressively enlarged 
systems and then estimating the energy of the corresponding infinite system by means of 
extrapolation to N — > oo. Certainly, approach (ii) results computationally more demanding 
than (i) but also more accurate, and we find a pressure difference of ~ 5 GPa between 
both resulting EOSs at the smallest studied volume (V = 2.50 cm 3 /mol). Nevertheless, 
this discrepancy by itself does not explain the large disagreement between our results and 
high-pressure data. In consequence, we turn our main concern to the characterization of the 
interatomic interactions. 

It is well-known that the structural and electronic properties of the atomic and molecular 
systems may experience important arrangements by effect of pressure.— i2£4ML Upon com- 
pression, overlappings between the electronic clouds within the system are promoted hence 
further correlations among the atoms (angular forces) emerge so as to lower their energies. 
In the case of solid 4 He, it has been suggested and tested within the Self Consistent Phonon 



formalism that 
sity.-^ In Ref. 
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;hree-body exchange interactions become significant with increasing den- 
Chang and Boninsegni included three-body effects into their high-pressure 
PIMC calculations performed with pair potentials, by computing the energy of several three- 
body interaction models over sets of configurations generated in their simulations (that is, 
perturbatively). In doing this, their agreement with experiments did not improve quanti- 
tatively, thus they suggested that higher order many-body contributions to the energy had 
to be considered. More recently, Herrero^ has adopted a similar but computationally more 
demanding approach to that of Chang and Boninsegni in which three-body interactions are 
explicitly included into the model Hamiltonian. For the case of a rescaled Bruch-McGee 
three-body interaction, Herrero's work provides very notable agreement with experiments 
up to pressures of ~ 52 GPa and at room temperature. 

In the present work, we propose a perturbative approach for quoting the many-body in- 
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teractions happening within highly compressed solid 4 He which are not accounted for by any 
atomic pair potential, and without increasing the computational cost significantly. Essen- 
tially, this consists in performing ab initio density functional calculations over sets of atomic 
configurations independently drawn from DMC simulations; subsequently, the energies pre- 
viously computed with DMC are corrected according to the average difference between the 
ab initio interaction energy of the all-electron-ion system and the pair-potential energy. In 
this way, many-body interactions of order two and higher are included perturbatively into 
the EOS without requiring from the knowledge of any additional two-, three-, four-, and so 
on, body interatomic potential. We show that proceeding in this manner the agreement with 
high pressure experimental data is at the ~ 5 — 10% level with relatively few computational 
extra cost. Truly, the approach that here we present for helium can be extended to the 
study of other quantum crystals upon high pressure for which accurate pair potentials are 
available. 

The remainder of this paper is as follows. Section [Til describes the methods that we have 
employed, the treatment of finite size effects and gives the technical details. In Section IHIl 
we present our results for the ground state of solid He at low and high pressure and yield 
further comparison with first-principles based calculations. Finally, in Section [IV] we analyze 
the pros and cons of the proposed perturbative scheme and conclude with the final remarks. 

II. APPROACH AND METHODS 
A. Diffusion Monte Carlo 

DMC is a ground-state method which provides the exact energy within statistical errors 
of many-boson interacting systems of interest.— >ili21 This technique is based on a short- 
time approximation for the Green's function corresponding to the imaginary time-dependent 
Schrodinger equation, which is solved up to a certain order of accuracy within an infinitesimal 
interval At. Despite this method is algorithmically simpler than domain Green's function 
Monte Carlo,— ^ it presents some (Ar) n bias coming from the factorization of the imaginary 
time propagator e~^~ H . Our implementation of DMC is quadratic,— hence the control of 
the time-step bias is efficiently controlled since the required At — > extrapolation is nearly 
eliminated by choosing a sufficiently small time step. The Hamiltonian, H, describing our 
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FIG. 1: Radial pair distribution function of solid 4 He at different molar volumes as computed with 
DMC. Curves are terminated at half the length of the simulation box (containing 180 particles in 
each case). 
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where mn e is the mass of a 4 He atom, the distance between atoms composing an i,j 
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the HFD-B(HE) Aziz interaction.— The corresponding Schrodinger 



equation in imaginary time {it = r), 
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(2) 



with E an arbitrary constant, can be formally solved by expanding the solution \P(R, r) 
in the basis set of the energy eigenfunctions {0 n }. It turns out that \&(R, r) tends to the 
ground state wave function <f) of the system for an infinite imaginary time as well as the 
expected value of the Hamiltonian tends to the ground-state value Eq. The hermiticity of 
the Hamiltonian guarantees the equality 

Vo|H|0 o ) (0o|H|Vr> 



E 



(0o|0o) (<ft#r) 



(3) 
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where ipr is a convenient trial wave function which depends on the atomic coordinates of 
the system R = {ri, r 2 , r^t}. Consequently, the ground-state energy of the system can 
be computed by calculating the integral 

E DMC = lim / E L (R) / (R, r) dR = E , (4) 

where / (R, r) = ^ (R, r) tpT (R) (assuming it is normalized), and El (R) is the local energy 
defined as -Ex(R) = Hip T (R) /tp T (R). The introduction of ipx (R) in / (R> T ) is known as 
importance sampling and it certainly improves the way in which integral (j4j) is computed 
(for instance, by imposing ipT (R) = when is smaller than the core distance of the 
interatomic interaction). 

In this work, the trial wave function adopted for importance sampling corresponds to the 
extensively tested Nosanow-Jastrow model^ 1 ^ 1 ^ 

N N 

V^j(ri,r a , ...jrrsr) = JJf 2 (r y ) JJgi(|r|-R||) , (5) 

i^j i=i 

with f 2 (r) = e~2(r) an d gi(r) = e~2 ar2 where a and b are variational parameters. This 
model is composed of two-body correlation functions f 2 (^) accounting for the two-body 
correlations induced by V^(r), and one-body functions gi(r) which localize each particle 
around a site of the equilibrium lattice of the crystal as given by the set of vectors {Ri}. 
Initially, the parameters contained in -^nj are optimized by means of variational Monte Carlo 
at some molar volume near equilibrium; however, as we have explored the system over a 
wide range of volumes we have repeated this procedure at some other selected points along 
the EOS. For instance, the optimized value of the parameters a and b at the molar volume 
20.48 cm 3 is 1.12 and 0.87 A~ 2 , respectively, while at 4.02 cm 3 results 1.15 and 3.06 A~ 2 . 
The parameters of the simulations, namely, the number of particles, N, critical population 
of walkers, n w , and time step, At, have been adjusted to eliminate any residual bias coming 
from them; their respective values are 180, 400 and 2.7- 10~ 4 K _1 . The parameter At has 
been reduced progressively with increasing density in order to provide numerical stability. 

B. Finite size corrections 

The description of an infinite system of interacting particles is obtained in practice 
through the simulation of a finite number of particles enclosed within a box. The difference 
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between the scale of the real and simulated systems can be overcome by enlarging the size 
of the simulated system so much as possible and applying periodic boundary conditions to 
it.— Even so, several corrections to the energies quoted directly from the simulation must be 
done if correlations of longer range are present. Certainly, these corrections arise from the 
fact that the maximum distance involving correlations in the simulation coincides with the 
length-scale of the particle container. The expressions for the potential and kinetic energy 
corrections AV tad and AT tad , assuming a certain cut-off length K max for the computation 
of the correlations (generally chosen as half the length of the simulation box), are : 

POO 

Av taii = 2nNp / g( r )V 2 (r)r 2 dr (6) 

" Umax 
POO 

AT tail = -AnNDp / g{r)V 2 \ni 2 (r)r 2 dr , (7) 

" Fornax 

where N, D = h 2 /2m and p are the number, diffusion constant and density of particles, 
and g(r), V 2 {r) and f2(r) the radial pair distribution function, pair potential and two-body 
correlation function entering the trial wave function, respectively. In the case of liquids, g(r) 
can be well-approximated to unity in equations and ([7]), and consequently, AV tad and 
AT tad turn out to be analytically accessible (standard tail correction -STC-). Nevertheless, 
in the case of solids such approximation could result rather inaccurate owing that the pattern 
of the radial distribution function is still oscillating beyond the cut-off distance (see Fig. [T]). 
In view of these facts and in order for the attained description of solid 4 He to be as precise 
as possible, we have estimated AV tad and AT tad also by means of VMC (variational tail 
correction -VTC-) through the relation 

AE tad = AT tad + AV tad = E™ MC - E^ MC , (8) 

where the superscripts in the energies refer to number of particles, N is the number of par- 
ticles used in the DMC simulations and E VMC = ( , 0t|H| , 0t) / (V'tIV't)- The limit iV — >• oo, 
equivalent to R ma x — > oo in equations ([H]) and ([7]), is reached through successive enlargements 
of the simulation box at fixed density (up to 900 particles) and further linear extrapolation 
to infinite volume. Indeed, this procedure results computationally affordable within VMC 
but not within DMC. In Fig. [21 we shown the asymptotic agreement between standard and 
variational energy tail corrections for infinite solid 4 He (1/N —> 0) within VMC. 
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FIG. 2: Variational energy per particle in solid 4 He at V = 21.35 cm 3 /mol as a function of 1/N. 
Filled circles correspond to total energy assuming STC energy tail corrections while the empty 
ones correspond to the total energy deduced directly from the simulation; both respective linear 
fits are coincident in the limit N — > oo 

C. Ab initio calculations and perturbative approach 

Density Functional Theory (DFT) is a first-principles quantum approach which has al- 
lowed for accurate and reliable knowledge of a great deal of materials and systems with 
exceptional computational affordability. A comprehensive description of DFT methods as 
applied to the modeling of condensed matter is given in recent books and reviews.— 1 ^ In 
DFT, the ab initio free energy of an atomic system, given the positions and charges of its 
nuclei, is expressed as a functional of the electronic density, n(r), as follows: 

E[n(r)} = T[n(r)} + \ ] J ^^drdr> + f>/ J^T]^ + 

N 

W)] + E i B 1 4 r ^ 
7?j 1 Rl " Rj 1 

where T[n(r)] is the electronic kinetic energy, Zj and Ri the atomic number and position 
of atom /, respectively, and E xc [n(r)] the electronic exchange-correlation energy (we have 



imposed l/Aneo and e = 1). The other terms in Eq.([9]) account for the Coulomb inter- 
actions between electrons, electrons and nuclei and nuclei. The Hohenberg-Kohn theorem 
states that the density n (r) which minimizes the functional J5[n(r)] corresponds to the true 
ground-state density of the system (thus £o({R}) = E[n (r)]) and that this optimal solution 
is unique. It is demonstrated that DFT is an exact electronic ground-state method, whereas 
the electronic exchange-correlation functional is not known for most of the systems. Conse- 
quently, some approximations for it must be introduced in the calculations. The most widely 
used models for E xc are the local density approximation (LDA) and the generalized gradient 
approximation (GGA), which have been parameterized by different groups. In this work, 
we use both Ceperley-Alder (CA) version of LDA— and Perdew-Wang (PW91) of GGA,— 
since a priori one can not discern confidently which is going to result more reliable for the 
study. A completely independent issue from the choice of E xc is the implementation of DFT 
that is used. This mainly concerns the way in which the electron orbitals are represented. 
Here, we use the projector augmented wave (PAW) framework developed by Bloch 5 - and as 
implemented in the VASP program.— 

The perturbative approach that we propose for correcting the DMC energies obtained 
with the pair potential V AzizI1 , consists in averaging the quantity AE = E ({R}) — 
Yli<j V AztzI1 (Rij) over sets of configurations drawn independently from the DMC simu- 
lations. According to this, the corrected energies result 

E'dmc = Edmc + (AE) DMC . (10) 

The many-body correction (AE) dmc includes two-, three- and so on many-body contribu- 
tions to the total energy as can be seen by invoking a many-body expansion of the ab initio 
ground state energy 

N N 

AE = E Q ({R})-J2V2 AzizII (Rij) = Y<V2(Ri j )-V 2 AzlzII (R lJ )+ V 3 (R ij ,R ik ,R jk )+--- . 

i<j i<j i<j<k 

(11) 

We note that the family of vectors {R} here refer to the positions of the atoms (nuclei) and 
not to the sites of the perfect crystalline lattice. It turns out that all the many-body terms 
composing AE are evaluated for any arrangament of the atoms as generated according to 
the Hamiltonian in Eq.(TjQ), and included into the total energy in a perturbative manner. 
Certainly, our many-body approach is not exact; firstly, it is noted that the full quantum 
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Hamiltonian of the system expressed within the Born-Oppenheimer approximation (FQ- 
BO) might be written down as Rfq-bo = T ion + E , where T ion corresponds to the kinetic 
energy of the nuclei, and that Eq.flT]) results a simplification of it, needless to be said, 
extremely accurate at low and moderate pressures. Nevertheless, using the DMC method 
for solving the ground-state of Hfq-bo remains a future goal because of the numerous 
intricacies encountered in the treatment of the electronic degrees of freedom (i.e. choice of 
the trial wavefunction and sign problem) and large computational cost involved. Therefore, 
instead of abording the full quantum problem straightaway, we have opted for a simplified 
but affordable strategy: add and substract V^ zlzI1 into Hfq-bo, solve exactly the part of 
the Hamiltonian embodying most of the two-body interactions and account for the rest by 
means of first-order perturbation theory. 

The ab initio calculations required for the computation of (AE)dmc have been performed 
on supercells containing 180 particles and with 2x2x2 Monkhorst-Pack k-sampling of the 
Brillouin zone^ and cut-off energy 478.0 eV; these settings ensure energy convergence to 
better than 0.5 K per atom. On the other side, our criterion for the convergence of correction 
(AE) dmc relies on the measure of its fluctuation, 5AE 2 = (AE — (AE)) 2 , over the same col- 
lection of DMC configurations than used for the average (-)dmc- Given a molar volume, we 
have requested y/ (5AE 2 ) DMC to be less than 1 K per atom, that is approximately 0.1% the 
total DMC energy obtained at large densities. In the process of drawing atomic configura- 
tions from the DMC runs, we have imposed the only constraint \El ({R-}) — < §|(-Ex)|, 
where El ({R}) is the local energy of the considered configuration and (El) the mean energy 
calculated over the population of walkers to which it corresponds. We have proceeded so 
for avoiding spurious configurations on the averages which otherwise are rejected within few 
steps in the DMC sampling. The number of atomic configurations required for the conver- 
gence of (AE)dmc has proved smaller than initially expected in all the studied cases: about 
15-25 were enough. This rapid convergence of the fluctuations a/ (5AE 2 ) dmc reveals that 
despite two-body interactions by themselves are not sufficient to attain reliable description 
of very dense solid 4 He they are still of leading relevance on it. 
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III. RESULTS 



A. Low pressure regime 

The EOSs of solid 4 He has been obtained by fitting a fourth-order polynomial to the 
DMC energies and subsequently performing the derivation respect to volume, 

where Vq is directly the equilibrium volume of the system and a, b constants. In Fig. [31 
we compare the DMC energies at volumes close to equilibrium (P(V) ~ 0), obtained with 
both STC and VTC, with experimental data; 1 As one observes there, excellent agreement 
between measurements and VTC results is provided, however, differences with respect to 
STC results are quoted in an almost constant upwards shift of ~ 0.30 K at positive pressure. 
Despite these discrepancies will practically vanish when expliciting both VTC and STC 
EOSs because of the energy derivative involved (as it will be shown shortly), it is noted 
that for other magnitudes which explicitly depend on the internal energy, as for example 
the enthalpy or freezing and melting densities, VTC and STC lead to different results. In 
the same figure, we also display previous theoretical calculations performed with Green's 
function Monte Carlo (GFMC) and the HFD-HE2 Aziz pair potential. 30 The GFMC points 
perfectly coincide with our results obtained with VTC, however, they disagree with the 
STC ones in the same manner experimental points do. Since GFMC and DMC are exact 
ground-state methods, energy differences between both approaches should be due only to 
the model interaction. Assuming that the treatment of finite size effects adopted in Ref. |30 



corresponds to the commonly used STC one, we may just conclude that both HFD-HE2 
and HFD-B(HE) interatomic potentials are likely to produce equivalent P — V curves at 
low pressures (likewise VTC and STC lead to practically identical EOSs) but not so total 
energies and other directly related properties. 

Fig. H] reports our results of the EOS of solid 4 He at T = in the volume range 21.0 < 
V < 7.5 cm 3 /mol (0 < P < 1.0 GPa). Curves obtained with VTC and STC are coincident 
as we anticipated in the previous paragraph. The parameters of the fit ffl2|) also displayed 
in Fig. H] (we provide the one obtained with VTC) are a = 9.11(6) K, b = 16.93(15) K 
and Vq = 25.04(4) cm 3 /mol (uncertainties are shown within parentheses). A glance at the 
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FIG. 3: Total energy per particle of solid He at low pressures as function of density (expressed 
in units of a = 2.556 A) computed with DMC and the HFD-B(HE) Aziz interaction. Results 
are obtained with VTC (•) and STC (A) and compared to experimental data of Ref. 1 (a) and 
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previous GFMC calculations (v) performed with the HFD-HE2 Aziz potential found in Ref. 
Error bars are smaller than the depicted symbols. 



plot reveals an excellent agreement between our results and experiments at low pressures, 
however, discrepancies become progressively larger as we move towards volumes smaller 
than 8.5 cm 3 /mol (P ~ 0.65 GPa). (For instance, at V = 7.76 cm 3 /mol our prediction 
of pressure overestimates the experimental value within ~ 10%.) It is worth noticing that 
the worsening of our curve roughly coincides with the interval in which the potential energy 
of the system becomes positive (see Table I). This fact indicates that the repulsive part of 
the HFD-B(HE) potential is probably too stiff. In the next subsection we will extensively 
deal with the shortcomings derived from the adopted model interaction, however now we 
continue with the description of other atomic magnitudes of interest that we have obtained 
at low pressures. 

The zero-temperature atomic kinetic energy of solid 4 He, E^, is an important (and chal- 
lenging) quantity to measure and compute since it evidences the singular quantum nature 
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FIG. 4: EOS of solid 4 He as computed with DMC and HFD-HE2 interatomic potential in the 
pressure range < P < 1 GPa. VTC and STC lead to identical curves within the statistical 
uncertainty and experimental data of Ref. 1 (points) is provided for comparison. 

of this crystal. It is well-known that the zero-point energy of solid helium is comparable in 
magnitude to its potential energy (cohesive energy), E p , and that the ratio between these 
two energies gives a qualitative idea about the relevance of anharmonic effects in the system 
(the larger E^jE v is, the larger anharmonic contributions would result).- 1 ^ From a compu- 
tational side, exact estimation of the expected ground state values of operators which do not 
commute with the Hamiltonian, as for instance the potential and kinetic energy operators, 
may be provided within the DMC scenario by means of the pure estimator technique.— 
In practice, this technique involves the introduction of additional weight factors into the 
customary DMC sampling which retain memory of the configurational replication processes 
occurring along the simulation. In order for our evaluation of the zero-temperature kinetic 
energy of solid 4 He to be as reliable as possible, we first determine the exact potential energy 
of the system by means of the pure estimator method and then we subtract it to the total 
energy (we note that within DMC the estimator of the kinetic energy is slightly biased by 
the choice of the trial wavefunction). In Fig. [5] we display our results for Ek and compare 
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them to low temperature data provided by several authors:— i22i2L£2i the overall agreement 
between them is excellent. In particular, we note the perfect agreement of our calculated 
value Ek = 24.24 (5) K at V = 20.87 cm 3 /mol with the very recent neutron scattering mea- 
surement of Diallo et a/.,— E^ xpt = 24.25 (30) K, performed at the same volume. In Table 
I, we enclose DMC results for the total, kinetic and potential energies of solid 4 He including 
STC at some selected volumes within the interval 22.60-8.0 cm 3 /mol. In Fig. [51 however, 
we have not refrained from including a further point at volume 6.70 cm 3 /mol for which we 
have shown that the description of the system attained with the model interaction seems to 
be not fully reliable. It should be mentioned that the treatment of finite size effects adopted 
in the calculations has little effect on the final values of Ek since the largest contribution to 
the total energy tail correction stems from the interatomic interactions. 

Another quantity of interest in the study of quantum solids is the atomic mean squared 
displacement, (u 2 ) , which is directly measured in x-ray diffraction experiments. In con- 
nection to this magnitude, the Lindemann ratio is defined as 7 = a/ (u 2 )/d, where d is the 
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FIG. 6: Left: Lindemann ratio of solid 4 He as function of volume. Our results are • and line 
0.058 + 0.0097 V as guide-to-eye, and experimental data of Ref. 



63 (A), 



6J (A) and |65j (V) are 



shown for comparison. Right: Debye Temperature of solid 4 He at T = as function of volume 



Our results are • and line 
comparison. 



and experimental data of Ref. 



6J (A) and 



63j (A) are shown for 



distance between nearest neighbors in the perfect crystalline lattice. As we pointed out in 
the Introduction, the zero-temperature Lindemann ratio of solid 4 He is uncommonly large 
(even if compared to other distinguished quantum solids like for example H2 which possesses 
7 ~ 0.18) as consequence of its light atomic mass and weak interatomic interaction. Using 
the pure estimator technique, we have studied the dependence of 7 with volume over the 
range 22.6-8.0 cm 3 /mol. We have depicted our results for 7 in Fig. [6] and compared them 
to experimental data of different authors, and again the overall agreement between them is 
remarkable. Once (u 2 ) is known, the Debye temperature of the system at T = 0, 0£>, is 
deduced straightforwardly through the relation Od = 9/i 2 /4mH e (u 2 ) . We have fitted our 
results for 0£> with the relation 

Q D = exp (^CiX* j . (13) 



,i=0 



where x = In (V/Vd) and which has been used previously to reproduce the density depen- 
dence of the phonon frequencies in solid H 2 and 4 He as well.— ^ Our optimal coefficients for 
expression (JEJ} plotted in Fig.[6]are: V D = 22.6166 cm 3 /mol, c = 3.21655 , c x = -2.23859 , 
c 2 = 0.122057 and c 3 = 0.319911 . (Additional point at V = 6.70 cm 3 /mol in Fig. El has 
not been used in the fit.) 
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y(cm 3 /mol) 


E DMC /N{K) 


E k /N(K) 


E p /N(K) 


22.60 


-6.51(2) 


21.36(6) 


-27.87(6) 


20.95 


-6.22(2) 


24.20(6) 


-30.42(6) 


19.34 


-5.50(2) 


27.63(6) 


-33.13(6) 


17.96 


-4.32(2) 


32.01(6) 


-36.33(6) 


16.76 


-2.50(2) 


35.24(6) 


-37.74(6) 


15.24 


1.63(2) 


42.63(6) 


-41.00(6) 


14.37 


5.25(3) 


47.09(8) 


-41.84(7) 


13.41 


11.11(5) 


53.66(9) 


-42.55(8) 


10.06 


68.80(5) 


89.90(9) 


-21.10(8) 


8.04 


192.45(5) 


133.00(9) 


59.45(8) 



TABLE I: Total, kinetic and potential energies per particle of solid 4 He including STC [Ebmc^ 
E k and E p , respectively) as computed with DMC and the HFD-B(HE) Aziz potential. Figures 
within parentheses account for the statistical errors. 

B. High pressure regime 

As we have already illustrated in Sec lHI Al the pair potential HFD-B(HE) performs excel- 
lently in the description of solid 4 He up to volumes of 8.5 cm 3 /mol, however, it monotonically 
fails to reproduce its EOS as the density is increased beyond this point. In Fig. [71 we ex- 
plicit the differences between measurements of Ref. 1 and [2] and our calculations performed 
with STC and VTC for finite size effects and over the pressure range — 57 GPa. The 
pressure difference between the EOSs obtained with VTC and STC at the highest studied 
density amounts to ~ 5 GPa, however this quantity results very small when compared to 



the discrepancy of both with experiments which is ~ 40% of the experimental value. T 
discrepancy is in overall agreement with the microscopic calculations of Refs. 



his 



37 and 



36 



Our results and other found in the literature,—* 3 ^ pose the need for considering higher 
order many-body effects present on dense 4 He instead of going in the search of improved 
pair potential models. As we pointed out in the Introduction, many authors have made 
efforts for elucidating the relevance of three-body and higher order (up to six-body) effects 
on the EOS of solid 4 He with assorted degree of accuracy and success.— MM. In this section, 
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FIG. 7: Equations of state of solid 4 He over the — 57 GPa pressure range as computed with DMC 
and HFD-B(HE) interaction. Experimental data of Ref. l] (•) 
comparison. 



and y (dashed line) are included for 



we present the P — V curves calculated within our proposed scheme for correcting the DMC 
energies obtained with pair potentials as described in Sec. Ill CI Just as we have explained 
there, all the many-body interactions not accounted for by V^ ztzI1 are computed with ab 
initio DFT and included into the total energy in a perturbative way, without requiring from 
the knowledge of additional two-, three- and/or higher order many-body interaction models. 

The fits to our results displayed in the P — V figures of this and next subsections have 
been performed with the Vinet relation^ 

2 r 

3 



P{V) = 3B 1 




cxp 



B n 



11 



V 
Vo 



(14) 



where Vo, -Bo and B' are the equilibrium volume, equilibrium isothermal bulk modulus 
{B = —VdP/dV) and equilibrium dB/dP, respectively. The experimental values of these 
parameters as provided by Ref. (2) are Vq Xp1 = 13.72 cm 3 /mol, B^ xpt = 0.225 GPa and 



B, 



'expt 



7.35 . We have also enclosed data points of Ref. LL in the plots for additional 



comparison with our estimations. The improvement of our EOS when considering many- 
body effects computed with the proposed perturbative approach is substantial (see Fig. [Sj). 
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DMC - HFD-B(HE) 
DMC + LDA correction 
DMC + GGA correction 
expt, Driessen 
expt, Loubeyre 




V (cm7mol) 

FIG. 8: Zero-temperature EOS of solid 4 He as computed with DMC and HFD-B(HE) pair potential 
and considering perturbative many-body corrections to the energy (solid and dashed-dotted lines 
mean LDA and GGA corrections, respectively). Experimental data of Ref. Q (•) and 2 (dashed 
line) are enclosed for comparison. 

For example, within LDA we obtain P = 54.83 GPa at volume 2.5 cm 3 /mol which is quite 
close to the experimental value 56.94 GPa and far below the non-corrected DMC result 



of 83.60 GPa. The parameters of the fit corresponding to this case are V c 



LDA 







7.77 



cm 3 /mol, Bq DA = 1.884 GPa and B' LDA = 6.66 , which as it is observed in Fig. [8] leads 
to constant underestimation of pressure within few Gpa respect to experimental data along 
the whole depicted range. On the contrary, the EOS obtained with GGA provides a notable 
description of the system near the equilibrium and few GPa above the experimental values 
at high pressures. Putting this into figures, V GGA = 12.93 cm 3 /mol, B GGA = 0.510 GPa and 



Bq* ga = 6.53 , which in fact result closer to the experimental values of Ref. 12| than the LDA 
ones. It is worth mentioning that the observed tendency of GGA (LDA) for overestimating 
(underestimating) pressure in our results is a well-known outcome in the field of ab initio 
simulations. 

Table II yields the values of the DMC energy and perturbative corrections for solid 4 He at 
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1/ (cnr/mol) 


E DMC /N{K) 


/ A 7~i\ TjDA /7\7"/TV\ 

( ae )dmc/ n ( k ) 


/A 7— i\ CI CI A /]\r/T/'\ 

( Ae )dmc/ n ( k ) 


10.06 


68.80(5) 


0.00(75) 


0.00(14) 


6.70 


404.55(5) 


-352.55(88) 


72.06(33) 


5.03 


1163.54(8) 


-813.08(55) 


200.43(39) 


4.02 


2444.11(12) 


-1407.99(50) 


232.38(35) 


3.35 


4294.67(15) 


-2165.61(61) 


43.53(50) 


2.87 


6728.33(38) 


-3113.77(67) 


-389.66(51) 


2.51 


9742.06(49) 


-4263.34(98) 


-1055.16(98) 



TABLE II: Calculated DMC energies and corrections (AE 1 ) dmc P er particle for solid 4 He at some 
selected volumes. Within the parentheses are the statistical uncertainties, which in the case of the 
corrections correspond to \J (SAE 2 ) DMC /N (we note that (98) = ±0.98 and (5) = ±0.05). 

some selected volumes. Separately, we have shifted all the LDA and GGA corrections a same 
amount for providing null contributions at the largest enclosed volume so as to facilitate the 
comparison between them. Certainly, this can be done without any loss of generality since 
the zero of the LDA and GGA ab initio energies and that of the HFD-H(BE) interaction 
do not coincide and we are essentially interested in the pressure. Two main conclusions 
can be extracted from the values ( AE) dmc m Table II: (i) corrections performed with LDA 
decrease monotonically with compression, not so with GGA, and (ii) GGA corrections are 
smaller in absolute value than the LDA ones. Since the proposed approach for correcting 
the DMC energies obtained with pair potential is perturbative, the conclusion (ii) concedes 
more reliability to the results obtained with GGA than with LDA. Indeed, a conclusive 
answer about whether LDA figures are or not too large would be best provided by second 
order perturbative theory, however this is out from our scope. In the next subsection, we 
will comment again on this issue by supplying further comparison between results presented 
here and others obtained by means of ab initio procedures. 

One known weakness of DFT calculations is that the usual approximations for E xc may fail 
to capture the essence of the long-range forces present in the system.—^ In the case of rare 
gases, the van der Waals (vdW) energy, which physically accounts for Coulomb correlations 
between distant electrons, has notorious relevance in the cohesion of the system. With the 
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Expt. 


DMC 


LDA 


GGA 


LDA(vdW) 


GGA(vdW) 


Vo(cm 3 /mol) 


13.72 


20.16 


7.77 


12.93 


7.06 


11.35 


Bo(GPa) 


0.225 


0.018 


1.884 


0.510 


3.072 


0.901 


B'o 


7.35 


9.85 


6.66 


6.53 


6.19 


6.13 


Pmax (GPa) 


56.94 


83.60 


54.83 


63.79 


52.42 


62.00 



TABLE III: Parameters of the fits performed with relation (I14p for the resulting EOSs. The headers 
on this first row correspond to experimental values of Ref. 2, DMC calculations with pair potential 
HFD-B(HE), DMC calculations with many-body corrections as obtained with LDA, GGA, LDA 
plus vdW interaction and GGA plus vdW interaction, respectively. P ma x is the value of the pressure 
obtained at the smallest studied volume 2.5 cm 3 /mol. 

aim of estimating the effect of this shortcoming in our corrections, we have added an effective 
two-body term accounting for the vdW interactions to the ab initio energy Eq. This term 
is expressed as 

v vdW (R) = f(R)^ , (15) 

where C 6 = -10130.639 KA 6 and f(R) = exp (- (D/R - l) 2 ) for R < D but f(R) = 1 for 
R> D with D = 4.392944 A (that is, as given by the HFD-B(HE) interaction), and it has 
been evaluated over the same sets of atomic configurations than used for the computation 
of (AE) 

dmc- Following this receipt, the many-body correction devised for energies E^mc 
now can be redefined as AE vdW = E ({K}) + J2i<j Vvdw{Rij) - J2i<j Vz(Rij). In Fig. El 
we plot the curves obtained with the correction AE vd w, and Table III summarizes the 
parameters of all the fits that we have performed (for LDA and GGA corrections including 
and not including vdW contributions). On one hand, a glance at the figure reveals that 
considering vdW interactions as explained above has in general little effect on the results, 
just a slight and otherwise expected lowering of the P — V curves within few GPa over the 
whole depicted range. On the other hand, the equilibrium properties of the system, as given 
by the parameters of the fits, change appreciably (see figures enclosed in Table III). This 
result seems to corroborate the accepted assumption that the effect of long-range interactions 
in the EOS of rare-gas solids becomes less important with increasing pressure.— 

In this subsection, we have not attempted to enclose any result for fee 4 He in the plots 
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FIG. 9: EOS of solid 4 He as computed with DFT and corrected for the zero-point motion of 



the atoms with the Mie-Gruneisen model (GGA and LDA). Curves presented in Sec. piIBp and 
experimental data of Ref. Q and 2 are also included for comparison. 

and/or tables since experiments indicate that hep is the only stable phase of solid helium 
at high pressures (~ GPa) and low temperatures, apart from a small fee loop region around 
melting between 15 and 285 K. Reassuringly, previous work based on first-principles cal- 
culations agrees to regard hep as the most energetically favourable zero-temperature phase 
of 4 He upon pressures up to 160 GPa.— In spite of this, we have carried out a series of cal- 
culations in highly compressed fee 4 He in order to check the predictability of our approach. 
Essentially, our results show no appreciable energy differences between the two phases within 
the numerical uncertainty. This outcome, however, appears to be not surprising since short- 
range interactions in helium are of leading importance, and the first and second shells of 
nearest neighbours in the fee and hep phases peak at practically indentical distances given 
a same density. 
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C. Comparison with ab initio based calculations 



Within the DFT formalism, the zero-temperature energy of a solid is usually written as 
a sum of two different contributions 

E (V) = E perf (V) + E mb (V) , (16) 

where E per f{V) is the energy of the perfect crystal (atoms frozen on their sites) and 
E V ib(V) = E harm {y) + E anharm (V) accounts for the motion of the atoms and is expressed 
as a sum of harmonic and anharmonic terms. In practice, E per f is obtained with standard 
DFT calculations and it involves affordable computations performed within one unit cell of 
the perfect crystal (apart from the summations involving periodic boundary conditions). On 
the other side, the estimation of E^ requires from some knowledge of the phonon-related 
properties of the solid of interest. In the case of heavy-ion crystals, the quasiharmonic ap- 
proximation in combination with finite displacement methods have allowed for an accurate 
description of the phonon frequency spectra.— 1 ^ The basic strategy underlying these meth- 
ods consists in distorting the perfect crystal by displacing certain selected atoms slightly 
from their equilibrium positions and then evaluating the atomic forces arising on the system 
by means of the Hellman-Feynman theorem and DFT. This approach, however, fails to re- 
produce solid 4 He since it provides negative (imaginary) phonon frequencies associated to its 
experimental stable phases at intermediate pressures.— Truly, the relevance of anharmonic 
effects in solid He makes the computation of its vibrational properties a tedious and compli- 
cate task which requires from approaches going beyond the harmonic and/or quasiharmonic 
approximations. It should be noted that within DMC this difficulty is circumvented since 
the phononic nature of the studied system is inherently cast into the method, hence further 
partition of the energy into static and vibrational parts is not required. 

In order to contrast our results presented in Sec. IIII Bl with other obtained with ab ini- 
tio based methods, we have computed the EOS of solid 4 He through the relation ([TBI) by 
evaluating P per f with DFT and P^b with the Mie-Gruneisen models 

P mb (V) - -— - — , (17) 

where Od is the Debye temperature, 7g the Gruneisen parameter which we approximate 
as 7g = — d\nuj D /d\iaV (with ftw D = k B Q D ) and R the gas constant. Indeed, we have 
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used for @d{V) the experimental relation provided by Driessen et al— since, as we have 
noted previously, the estimation of this or any other vibrational property of solid 4 He at 
T = would result not trusty with customary ab initio strategies used in the study of 
normal (not quantum) solids. The resulting P^ increases monotonically with compression 
and for instance it represents about 15% of the total pressure of the system at volume 2.5 
cm 3 /mol, thus it must not be neglected. In Fig. [9j we show the EOSs obtained with the 
already explained procedure using both LDA and GGA approximations for the exchange- 
correlation energy; also we include the P — V curves quoted within DMC and corrected 
for the vdW energy and many-body interactions, and experimental data. As it is observed 
there, differences between both LDA and GGA perturbationally corrected curves and their 
ab initio counterparts are quite small. Moreover, these discrepancies are likely caused by 
the treatment of the long range interactions (vdW energy) described in Sec. IIII Bl and the 
approximation adopted for P^b- This result is stimulating since it demonstrates that with 
the approach presented in Sec. Ill Cl one can obtain accurate EOSs for dense solid helium, or 
equivalently for any other light quantum solid, in excellent agreement with those which would 
be obtained by means of first-principles approaches but with the benefit of not requiring 
from the computation of the phonon dispersion curves of the crystal or experimental data. 

Now we turn our attention to the concern posed over LDA in the previous subsection. As 
we noted there, a glance at Table II might lead to the conclusion that the LDA corrections are 
too large so as to be considered perturbative on top of the DMC energies (not so the GGA 
ones). Indeed, we do not dispose of a fair criterion for accepting or rejecting corrections 
in basis to their size, and in our opinion this is the most important shortcoming of our 
approach. Nevertheless, appealing to the good accordance between the LDA P — V curve 
corrected for the atomic zero-point motion and the DMC one corrected with LDA, we may 
feel quite confident about the reliability of the latter. 

IV. DISCUSSION AND CONCLUSIONS 

In the Introduction we pointed out that Diffusion Monte Carlo is among the best suited 
methods for studying quantum solids. In the case of bosonic systems, this method provides 
the exact ground state energy and related properties without dependence on the choice of 
the trial wave function, which otherwise is related to the computational efficiency Here, we 
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have proved the excellent performance of DMC with the Aziz pair potential HFD-B(HE) in 
characterizing solid 4 He at low pressures (P < 0.65 GPa), by estimating its EOS, kinetic 
energy per particle and Debye temperature at different volumes, and comparing our results 
with experiments. Especial attention has been paid to the extend of finite size effects in 
our results. To this regard, we conclude that customary strategies devised to correct such 
effects based on the approximation of the pair radial distribution function to unity beyond 
certain cut-off distance, may result accurate enough for the derivation of the EOS but not 
so for the assessment of other magnitudes like the energy. 

On the other side, solid helium under high compression (as most of the materials) un- 
dergoes important arrangements on electronic structure which lead to the appearance of 
angular correlations among the atoms. 41 This circumstance makes necessary to consider not 
only atomic pair interactions but also higher order many-body ones when investigating this 
crystal upon high pressure. Nevertheless, within DMC the minimal inclusion of three-body 
interactions on the model Hamiltonian has already the effect of drastically slow down the 
simulations. Furthermore, even in the supposed case computational cost was not a prob- 
lem, first we should know much better than now the analytical form of these many-body 
interactions or alternatively to be able to devise them (which actually may result puzzling). 
According to this occurrence, ab initio methods emerge among the best candidates for 
quoting such contributions since they do not rely on potential models and, in general, are 
computationally affordable. However, fully ab initio analysis of crystals requires from the 
knowledge of the phonon-related properties, and for the case of solid 4 He and other light 
quantum solids this is by no means a straightforward task. 

In this work, we have presented an approach for the study of dense solid 4 He at zero 
temperature which combines the versatility of the DMC method with the accuracy of ab 
initio calculations. On one hand, we naturally circumvent the calculation of the vibrational 
properties of helium thanks to the DMC strategy, and on the other, we account for the 
many-body interactions having place on the system by means of DFT. However, the way 
by which we enclose these many-body contributions to the DMC energy is not exact but 
perturbative and we do not dispose of rigorous tests for quoting the errors included on these 
corrections. Concerning results, we have yielded the EOSs corrected in this manner within 
LDA and GGA for the exchange-correlation energy and they have proved in fairly notable 
agreement with experiments over the pressure range — 57 GPa. Specifically, GGA provides 
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a better description of the crystal near equilibrium than LDA. Further comparison of these 
curves with EOSs obtained trough DFT and corrected for the atomic zero-point motion by 
means of an approximate model, comes to support the reliability of our approach. 

It should be mentioned that the zero-temperature scheme proposed in this work is also 
well suited for the study of other light quantum solids upon high pressure, like 3 He, H 2 , D 2 
and Li, for which accurate pair potentials are devised. Certainly, a further and promising 
improvement of the present framework would consist in going beyond the perturbative ap- 
proach. This could be achieved by proper coupling of the DMC and DFT methods, as for 
instance, by considering the ab initio energy of the system within the branching weight of 
the DMC algorithm. Work on this direction is in progress. 
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